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ANALYTIC SOLUTION OF THE SPENCER-LEWIS 


ANGULAR-SPATIAL MOMENTS EQUATIONS 


I. INTRODUCTION 


The major accomplishments of the project were 


1. A closed form solution for the angular~spatial moments of 


the Spencer~Lewis equation!/? 


2. A computer program that generates benchmark solutions for 


electron densities as a function of position and path length due 


to monoenergenic plane sources in infinite media. 


The original aim of the proposed work was to use benchmark 
calculations to aid in the development of numerical codes for 
It was proposed to 


determine energy deposition profiles using Spencer's method of 


solving electron transport problems. 


moments*, Such solutions would then be used to evaluate some 


of the numerical approximations used in the streaming ray*~* (SR) 


and SN®“* methods. Spencer's technique consists of reconstructins 


energy deposition profiles from angular-spatial~path length 


moments. 


During the course of this project it was discovered that 


the angular-spacial moments could be determined directly, 


thereby eliminating the need to take moments in path length. At 


this point it was decided to change the emphasis of the work. 


Instead of concentrating on the application of existing 


benchmark methods, the primary objective become the development 
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of new benchmark solutions based on the new closed form 


solution for the angular-spatial moments. 


The major effort was devoted to calculating the angular- 


spatial moments. The remaining work involved reconstructing 


electron density distributions from the moments. Several 


sample calculations have been carried out with monoenergetic 


plane sources of electrons in infinite media of aluminum and 


carbon. 


An important feature of the new solutions is that they 


provide the electron density distribution in position and path 


length rather than an integrated quantity such as the energy 


deposition profile. With a little more programming it appears 


that the angular electron density could be determined as well. 


qT. 


THEORY 


This section begins with an outline of Spencer's method of 


moments. It is then shown that the set of coupled ordinary 


differential equations for the angular-spatial moments can be 


solved directly, thereby eliminating the need for path length 


moments. Finally some of the standard techniques for 


reconstructing the electron density from its moments are 


discussed. 


Tl. A. Spencer's Method of Moments 


The electron distribution in space, angle and path length 


can be obtained by solving the Spencer-Lewis equation. In one 


» 
b 
PB 


dimension this equation is!/2 


(2 +4 2 +a) 9(XSoy) 


ce 


= 29 | o(y'on) o(X-Spu') du’ + Q(xX,S,y), (1) 


ba | 
where 

4(X,S,pu) = the density of electrons as a function of 

position x, path length s, and direction 4, 

0 = total interaction cross section, 

o(y'+y) = differential scattering cross section, 

Q(x,S,y) = the fixed source of electrons. 
We will assume that the medium is infinite and homogeneous, that 
initially there are no electrons present, 

o(x,O- 4) =O , (2) 
and that Q is a function of the form 

Ax,Su) = 6(x) 6(s) Ay) - (3) 
It is convenient to measure x and s in terms of the maximum 
electron range. Then the region of interest is defined by -l < 
x<lando<s<l. 


To solve Eq. (1), let 


0(%S0u) = x 262) yy (xes) Plu), (4) 
L= . 


o(Spy'ey) = > TOC (5) 
L= 


O(X Sry) = 6(x) é(s) 7 ith ae Po(u) ’ (6) : 
2=0 na 


where y, is the cosine of the angle between directions , and ,', 
and the P, are the Legendre polynomials. 

Using the addition theorem and recurrence relation for 
Legendre polynomials, and Eqs. (1) and (3) to (6), one obtains 


a 2 $94+1(X,S) = 6(x) 8(S) qg, 2 = 0, 1, 2,2... (7) 


Since the electrons originate at x = 0, it follows that 
o(X,Sep) = 0, |x] >s e (8) 
Thus multiplying Eq. (7) by x" and integrating over the interval 


-1 <x <1, yields the result 


[4 


[ds + a(s) - o4(5)]o ants) * aa oe-1,n-1(S) + 


(241) 6941,n-1(S) = 6n,9 6(S) Gp, Met 20, 1, 2,...-, (9) 


where 


%en(s) = x" 9,(x,s) dx . (10) 


The s dependence of the cross sections appearing in Eq. (9) 


is reasonably well modeled? by the relation, 


a: Dl 


=_ooo a 
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£ 
a(S) - oy(s) % 3 ' (11) 
where 
f, = 0(0) ~ a4(0) « (12) 


Using Eq. (11), multiplying Eq. (8) by (1-s)P*!, and 


integrating from 0~ to 1 yields Spencer’s result, 


perl pr | 
Ps én0% nee ns ey (2+) o i el 


ire Sadia iM | = 
in = Srltfy * (prl+f,) (2041) mee (y) 
where 

Pp 3 

* yn = | (1-s)P éen (s)ds . (14) 


Surface terms do not appear in Eq. (13) because the factor 
(1-g) (Pt) » Pp >-l vanishes at the upper limit, and as shown in 


Section IT. B., ¢gn(s) vanishes at the lower limit. The oP 
A 


. P 


>~1 are easily generated from Eq. (13). 


Il B. Closed Form Expression for the Angular-Spatial- Moments 

The simple recursion relationship for the vy Eq. (13) was 
made possible by assuming that the cross sections vary inversely 
with the residual electron range, Eq (11). It is shown here 
that this same assumption enables one to solve Eqs. (9) 
directly. 


For the case nm=0 we have 


i”, 


d fy 
(a5 + Tes) 4, ,(8) =q, é(s) . (15) 


From Eqs. (2), (4) and (10) it follows that 
o¢n(0-) = 0 fA n=0, l, 2, eee e (16) 
Thus, the solution of Eq. (15) is easily shown to be 
Cae a(s) (I-s)f, , 20 , 
9 =  ¢ (17) 
(_% e(s) , £20 , 
where @(s) is the step function. 


For n > 0, 


(f + ft) @en(S) = Spy | t62-1,n-1(8) + (£+1) sem (18) 
To solve this equation it is necessary to know the initial 
values of the $,7,(s) as s approaches zero from the right 
rather than from the left. 

Using Eq. (10), 


1 


lim ¢,n(s) = lim | x 9,(x,s)dx . (19) 
s+0* s+0 Ise 


From Eqs. (4) and (8) 


64(x8) =0, |x] >Ss, (20) 
therefore 
69n(07) =0Q, n>oO. (21) 


In spite of Eq. (20), 9,,(0*) # 0, as seen by Eq. (17). Thus, 
@,(x,0) must contain a delta function centered at x = 0 of 


strength qy. This should not be surprising in light of Eq. (6). 


After multiplying by (l-s)~£, Eq. (18) can be easily 


integrated to yield 


Lop—1 me1l(©) + (2+1)694+1,n-1(t) 


s) = n(l-s)f dt 22) 
Pant , (2g+1)(1-t) fy : 
i n>0O 
where use has been made of Eq. (21). ee 
Consider now the anzatz: 
n Yni ro 
oon (s) = $ > Agnizy (-s)(Eari + [il29) (23) : 
i=Ipn j=0 
where 
; 
Ign = max -n,-& ’ (24) 7 
((ne]il)/2, reli] even, 
Jni = { (25) 
{_ (rlil-1)72 , mil odd, 
and the Aynij are constants to be determined. ic 
Eqs. (17) and (23) are consistent provided 
Agoee = qa: (26) 


After some rather tedious algebra it can be shown that Eqs. (22) 


and (23) are equivalent provided the A,pij, n > 0 are given 


recursively by 


w. 


@I. 


Cynij + Egnij + Tyym2Sism2, jf, 
\ 


a eee Iun<i<tI 70 
snij an S tene2e Ws 
Agnij = 4 i . 
| Egnij PSPS oe TRO ; 
| Cenij + TinSis0 , JO , a 
n  Jni' 2 
FS aay seo . 
i'=Iy, j'=0 
io : 
where 
ngAy— i+] ,j- 
: tAreL meLitl inl , O<ic<m2, 340 , ° 
(2241) (£y4i+|i|+2j-£,) 
Cenij ib 
-__MAge1 nL itl, Ign Si < 0 
(2n+i)(Eppitlilt2j-fe)  ’ i 
and 
__£+1)Ag+1 ,n-1,i-1,5 l<icn 
9 
(224+1)(£ p44+|i|+2j-f2) se il 
eis : (29) 


i, 
n( £+1)Ag+1 rr] i-l, j~1 


( (22+1)(£ 945+ i|+2j-£2) 
TI. C. Computational Considerations 


ne Tree 1 C0" 290, 


Equations (23) to (29) constitute an exact solution for the 
o_en(S). We now discuss some of the details of their numerical 


implementation. 


Calculation of the Agni}: 


From Eqs. (23) to (25) it is evident that the A,ynij comprise : 


a very sparse array. On the other hand, as seen from Eqs. (27) 


to (29) the Agni; depend only on the A,,n-1,ij- Thus to avoid 
excessive storage requirements, the A, ,p-1,ij are discarded 

| after the Apgnij are determined. In addition, for a given n 

value the Aynjj are stored in a one-dimensional full array whose 


index, k(2,n,i,j) is given by 


g-1 n 
k=} YF pin) 
; g'= V=Ip'h 
| i-] | 
: rs S  pirtl) +541 (30) 
i'=Iyn 
It is possible to evaluate the summations analytically, 


f however the algebra is quite lengthy and will not be reproduced 


here. The final expression for k can be presented most 


compactly if the rules for integer arithmetic (any quotient is 


truncated to an integer value) are assumed to apply. In this 


case 


k=kjp+kg+j+1 (31) 


where 
a 


QO, 2£=0, 


n({n(122-2n-9)+36 2-10} 
ky =. ke RE, 7 mace e L>n, (32) 


i“ 6n2£(m+2+3)4+(2-1) [(2-1)(1-22)+26] 
24 


{ 0, vf = Iune 


(2((i-Tgn) (1 )=1]-(i-1) 12) 
4 , 
2[(i-L yy) * (ntl) +i] +(i-1)2-12, 


{ 4 
Mae 


kg =< 1+ Inn <i < 1, (33) 


> ae lee 


Ti. RECONSTRUCTION 

The electron density distribution can be obtained via Eq. 
(4) if the 9,(x,s) are known. At least on physical grounds, one 
would expect the @,(x,S) as functions x to belong to L*-1,1) 
space for any s value. Therefore, in theory it should be 
possible to reconstruct (with convergence in the mean) 9,(x,S) 
from its moments ,7(s), n = 0, 1, 2, ... . Since the 9,,(s) 
are determined recursively, in practice one must approximate 
$,(x,S) from a finite number, N of its moments. We consider 
three straight forward techniques. 
1. In the first method (and also in the second) we exploit the 


fact that there is a wave front located at |x|=s [see Eq. (20)] 


and put 
CS ams 
‘> SF tmls) PZ), [x] <8 
m=0 
o,(x,s) = < (34) 
0) ’ IxJ >s , 
ee 
where 


s 
Hym(s) = = | PCS) oy(xssldx , (35) 
-s 


and the Legendre polynomials, Pm(=) form an orthogonal basis on 


(-s,s]. If£ the coefficients, hm, are defined such that 


™m h (> 
Hem(s) = y = x" p(x,s)dx . (37) 


no § leg 


From Eqs. (10) and 


m 
eon 

Hem(s) = 2 <7 oan(s) « (38) 
rn=0 § 


2. The second procedure to be considered is mathematically 


equivalent to first provided all calculations are carried out 


with infinite precision. Therefore the numerical discrepancies 


between the two techniques should give some indication as to the 


severity of computer round off errors. 


In place of Eq. (34) we set 


am | 
> Gym x”, Ixl<s , 
m=0 
(XS) a (39) 


QO, [xj >s , 


— 


Multiplying Eq. (35) by x" and integrating from -s to s gives, in 


matrix form 

3, 7X6, (40) 
where G, and §, are N-dimensional column vectors whose m th sig 
components are Gym and 9,,,(s) respectively, and X is anN x N 


matrix whose n,m th element is given by 


s | a , Trtm even, 
ae xmim dx =~ gy ,ontm odd. os) 
Inverting Eq. (40), ‘ 
& =xls,. (42) 


3. The final reconstruction method is identical to the second 
except that the location of the wave front is not explicitly 
built into the procedure. The equations for this case are the 


same as those for method 2, except that Eq. (39) becomes 


N 
@_(X,S) = > Gem xm (43) 
m=0 


and Eq. (41) becomes 


Tm ’ nt+m even, 


Q , mm odd. (44) 


While methods 1 and 2 are theoretically identical for all 
values of N, this last method should agree (convergence in the 


mean) with the former two only in the limit Nee, 


= oy eee 


Iv. NUMERICAL RESULTS 


A computer program was written which determines the 
moments oF , and the total electron density, (x,s) defined by 


o(X,S) = 29 o(X,Sep) Gy = 9,(X,S), (45) 


as 
where the second equality follows from Eq. (4). The oa are 
determined both by Spencer's method, Eq. (13), and also from Eq. 
(23). The total density distribution is reconstructed by all 
three methods described in Section I. C., namely; Eqs. (34), 
(39), and (43). 


Iv. A. Cross Sections 

The moments method appears to be restricted to infinite 
media problems, and to cross sections of the form specified in 
Eq. (11), or the somewhat more complicated expression?,?® 


af, 


o(S)~a,(s) = Tics)(i-s-a)’ (46) 
where a is an adjustable parameter. Due to these limitations, 
moments calculations are probably more useful today for 
benchmarking more flexible numerical techniques, such as Monte 
Carlo'* or discrete ordinates*“*, than for producing realistic 
electron transport results. For this reason our results are 
based on Eq. (11) rather than the more complicated Eq. (46). In 
References 2 and 9 Spencer gives the counterpart of Eq. (13) for 


the case when the cross sections satisfy Eq. (46). 
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Although Eq. (11) is of limited validity it does permit the 
cross section for the incident (s=0) electron energy to be 
specified to any desired accuracy. The model assumed in all -- 
calculations is the following expression for screened Rutherford 


elastic scattering?, 


2 
= 2x2(Z+1) Etl | | 1 
o(s=0, Uy) — Na r,? Ee E,t2) Ti¥n=,)? e (47) 


where 
E, = kinetic energy of source electrons (s=0) in mc2 units, 
Z = atomic number of the transport medium, 
A #= atomic weight of the transport medium, 
Na = Avogadro's number, 
YX, = e2/mc? (classical electron radius), 
ue = cosine of the scattering angle, 


and the screening constant, n, is given by the Moliere formula’ 


ro2v3 Tro. 7 z\) (Eth)? 1)? | 
n= 05 | “aS5(137) [exer ina +206 by Bera | “4° 
j 


As specified in Eq. (48), the units of the cross section are : 
(cm?/g). To be consistent with our former development where 
distances are measured in range units we let 

o(s=0, y,) + o(Ss0, u,) * S(E,) , (49) 


where S(E,) is the maximum range in (gm/cm2) of electrons of 


energy E,. 


Iv. B. Moments 


According to Eq. (23) and the definition of Pn, Eq. (14) 


n 
= 2 a (50) 


islyn j=0 fei +|il+2}+p+1 


if £24; +/i]+2j+p+l > 0 for all 2, i, j, p and diverges otherwise. 
For the case 2=i=j=0 the denominator inside the double summation 
becomes p+l (f= g-0, = 0) and therefore, Eq. (50) has the 
restriction p >~1, which is identical to the restriction placed 
on Eq. (13). The oP, were calculated for :,p = 0, 1, 2, 3 and n 
= 0, 1, ... 40 with an isotropic (q, = 6,,) source of electrons 
and initial energies, E, = 1.0 MeV and 0.1 MeV. As_ shown in 
the left half of Tables I and IT the results obtained from &qs. 
(13) and (50) are in agreement to at least five significant 
figures provided that the calculations for Eq. (50) are carried 
out in double percision (on a CYBER 175 computer). However, Eq. 
(50) does appear to be sensitive to round off errors as seen in 
the second column of results in Tables I and II, which were 
obtained using single percision. 

An outside verification of the new formulation would be 
desirable. Unfortunately, Spencer's published moments values” 
are based on Eq. (46) rather than Eq. (11). Nevertheless, as a+e 
the two cross section expressions become identical. Comparison 
of moments determined from Eq. (50) and Spencer's published 


results are shown in Table III for an isotropic plane source, 14 
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=0, p= +. + owe i and three different values of «. The 


case a = 1.288 corresponds to 2.0 MeV electrons in aluminum, a = 
Gi 35.45 corresponds to 0.1 MeV electrons in aluminum, and a = 
107.8 corresponds to 0.025 MeV electrons in carbon. As 


expected, the agreement improves as |q| increases. 


Iv. C. Reconstructed Electron Densities 

The electron density »,(x,s) was reconstructed from the 
den(S), n = 0, 1, 2,2. N for 1.0 MeV and 0.1 MeV isotropic 
electron sources in aluminum. All three methods described in 
Section I. C. were used. For small N values the first two 
methods give better results since there was no need to 
represent the rapidly varying density, ¢,(x,s) across the wave 
front. On the other hand, with increasing N the third method 
appears superior at least for small s values. As s+0 the 
monomials x", n = 0, 1, 2, ... N begin “loosing” their linear 
independence on the interval [~s,s], that is, the matrix defined 
by Eq. (41) becomes more ill comditioned as s+0. Although a 
special matrix inversion routine developed for ill conditioned 
matrices was used, it eventually broke down for s small enough 
or N large enough. 

No difficulties were experienced inverting the matrix 
defined by Eq. (44). However, method 3 was subject to round off 
errors in the determination of the 6,,(s). The best results for 


the 1.0 MeV source electrons were obtained with N = 34 while 


the best results for the 0.1 MeV source occurred with N = 30. 
This conclusion is based on qualitative observations of the 
a results for larger N values, where the influence of numerical aa: 
errors is quite apparent. No independent solutions are 
presently available to make quantitative comparisons. It is 
i planned to repeat these calculations using the streaming ray 
numerical technique so that an independent verification will be 
possible. 
E The results of the electron density calculations are shown : 
in Tables IV-IX. The agreement between the three methods is 
quite good for large s although the breakdown of methods 1 and 
2 at small s values is apparent. This breakdown could be 
avoided by decreasing N with s. In fact, with the right choice 
of N, methods 1 and 2 should do at least as well as method 3 
since there is no need to numerically represent »,(x,s) = 0 in 
the region beyond the wave front. 

The electron densities determined by method 3 for plane 
isotropic sources in aluminum are plotted in Figs. 1 and 2 for 
E, = 0.1 MeV and in Figs. 3 and 4 for E, = 1.0 MeV. Figures 5 
and 6 represent the same physical situation as Figs. 3 and 4 
with the modification ¢(s=0) + 1.0 (range unit)“, Figures 1, 3 
and 5 give 4,(x,s) in the region .025 < s < .975, .0125 < x ¢ 
.988. The oscillations in the region {x| > s, s small, are due 


to the problem of trying to represent a delta function with a 
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finite number of polyncmials. This gives one incentive for using 


method 1 or 2. Alternatively, method 3 could be used to 

reconstruct the density of collided electrons which would be ete 
added to an analytic expression for the density of uncollided 

electrons. (The discontinuities at the wave front including the 

delta function at s=0 arise from uncollided electrons.) In 

order to cut out some of the nonphysical oscillation in the 

region |x| > s, s small, Figs. 2, 4 and 6 give 9,(x,s) for 0.25 < 

Ss < .98l. 

It is reassuring that method 3 is able to locate the wave 
front. This is more apparent from Figs. 5 and 6 then from Figs. 
1 to 4. Only uncollided electrons appear at the wave front. 
From Eq. (11), o(s) + # as s + 1. Thus the wave front 
disappears in Figs. 1 to 4 as s increases. To verify that 
method 3 correctly locates the wave front the mean-free-path 
of the electrons corresponding to Figs. 5 and 6 was set 


equal to 1.0 (range units). 


V. CONCLUSIONS 
A closed form solution for the angular-spatial moments of 
the Spencer Lewis Equation has been found. The validity of the 
solution was checked in two ways. First, the 9,7(s) were 
integrated to yield the oP. These agree to at least five 
significant figures with those obtained directly from Spencer's ; 


recursion relation, Eq. (13), and they are in good agreement with 
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Spencer's published results®’ for large a values (Eq. (50) 
corresponds to the case a = #). Second, the electron density, 
6,(X-S), was reconstructed from the %en(s)« Although no 
independent calculations are currently available for comparison, 
the results look reasonable and the location of the wave front 
is correctly predicted. 

Since the solution for the 9,,(s) are valid only for the 
approximate s dependence of the cross sections specified by Eq. 
(11), and for homogeneous infinite media with plane sources of 
electrons, their use for generating solutions to practical 
engineering problems is limited. On the other hand, they enable 
analytic benchmark solutions for @(x,s) and in theory for ¢(x,s,u). 
This is expected to be the main benefit of this project. 

Several interesting observations are possible as a result 
of the analytic expressions for the 4,,(s). The exponents in 
Eq. (23) are all positive except for the case 2=i=j=0 for which 
the exponent is es Therefore, 

\ 


Limit gen(s) =< 
sel 


Aonoos 2=0, 
0, otherwise , (51) 


implying [see Eqs. [4) and (10)] that the electron distribution 
becomes isotropic as s+l. The physical reason for this is that 
the differential scattering cross section blows up as sl 
{see Eq. (11)]. 

From Eqs. (17), (23) and the last line of Eq. (27), 


' Ap r=0 


Oo, m0 (52) 


limit 9,, (s) = ; 
in 
s+0F 
which is in agreement with Eq. (21). 
The total number of electrons in the medium for a civen s 


value is given by 


$,(X,S) dx = 9,,(s) = q, e(s) (53) 


i 
where use has been made of Eqs. (10) and (17). 

Therefore, 49(x,S) does not tend to zero at the electron 
range, that is, as s+l. In spite of the infinite scattering 
cross section, there is no mechanism in the model forelimi- 
nating electrons, and 4,,(S) remains constant. (For s > 1 the 
scattering cross section, Eq. (11), goes negative and &q. (22) 
predicts complex solutions.) 

As a final point, Figs. 1 and 2 for 0.1 MeV electrons are 
practically identical to Figs. 3 and 4 for 1.0 MeV electrons. 
Some of the similarity is due to scaling and the fact that 
distances are given in electron range units. Further studies 
are under way to determine if any additional explanations can be 


found. 


Iv. Suggested Follow-on Work 


From the results obtained so far, it appears that the 


following additional work should be carried out. . 


A. Modify the computer code such that the distributions of 
uncollided electrons is determined analytically, while moments 
1 reconstruction is used only for the density of collided ak 
electrons. This should speed convergence for Methods 1, 2 and 3 
of Section II.C and eliminate the oscillations in the region 


} |x|>s obtained with Method 3. 


B. Compare the electron densities determined with the modified 
code of Section IV.A with those obtained by an independent 2 


method such as the streaming ray numerical technique. 
C. Equations (23) to (29) constitute an exact solution of the 
> Spencer-Lewis angular-spatial-moments equations and should be 7 
reported in a scientific journal. This would require writing a 
F suitable introduction to Sections II to V which are already in = 
A journal format. In addition, some modifications based on the 
results of VI.A and VI.B would also be needed. 
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